R Notebook

library(here)
library(readxl)
#library(papeR)
#library(outliers)
library(kableExtra)
#library(DataExplorer)
library(lubridate)
library(forecast)
library(nlme)
#library(nortest)
library(ggfortify)
library(dygraphs)
#library(seasonal)
#library(seasonalview)

library(nonlinearTseries)
library(fNonlinear)
library(fGarch)
library(TSA)
library(tsDyn)

library(tidyverse)

Rutas

raw_data <- here("data", "raw")
interim_data <- here("data", "interim")
final_data <- here("data", "processed")

Leyendo base de datos

base <- read_excel(paste0(raw_data,"/Base Datos.xlsx"), 
    col_types = c("text", "numeric", "numeric"))

head(base)

Extrayendo la base de colones y de dolares

colones <- data_frame(date = base$`Activo neto`,
                  col =  as.double(base$CRC)) %>%
  mutate(date = ymd(paste0(date, "-01"))) %>%
  mutate(year = as.factor(year(date)),
         month = as.factor(month(date)))

head(colones)

Convirtiendo los datos a series de tiempo

# Colones

# 6 periodos para utilizar a modo de validación
colones_val <- colones %>% 
  slice_tail(n = 6)

# Serie completa para el analisis exploratorio
colones_full <- colones

# 10 años de datos para modelar
colones <- colones %>% 
  filter(date>"2007-12-01" & date<="2020-01-01")

# Objetos ts tanto para la serie completa como para la serie de modelado
colones_ts_full <- ts(colones_full$col, start = c(2001,2), frequency = 12)
colones_ts <- ts(colones$col, start = c(2008,1), frequency = 12)
colones_ts_val <- tail(colones_ts_full, 6)

1 Comprobacion de linealidad de la serie

1.1 Prueba de linealidad de la media

# Null hypothesis: Linearity in "mean"
tnnTest(colones_ts, lag = 1, title = NULL, description = NULL)

Title:
 Teraesvirta Neural Network Test

Test Results:
  PARAMETER:
    lag: 1
    m|df: 2
    t-lag-m|df: 142
  STATISTIC:
    Chi-squared: 7.9259
    F: 3.989
  P VALUE:
    Chi-squared: 0.01901 
    F: 0.02063 

Description:
 Thu Sep 10 20:06:28 2020 by user: 

La hipótesis nula del Teraesvirta Neural Network Test es que la media de la serie es lineal. Al rechazarse la hipótesis con una confianza del 95% se puede decir que la serie es no lineal.

1.2 Prueba para determinar si es caotica o no

options(max.print=1000000)
rqa.analysis=rqa(time.series = colones_ts, embedding.dim=1, time.lag=1,radius=0.01,lmin=2,vmin=2,do.plot=TRUE,distanceToBorder=2)

2. Inspección

# Una forma visual de empezar a revisar si existen o no clusters de volatilidad 

colones_ts_nd <-diff(log(colones_ts))


colones_ts_nd2<-colones_ts_nd-mean(colones_ts_nd) # Es el cambio relativo ajustado por la media en el tipo de cambio 
#colones_ts_nd2

colones_ts_nd3<-colones_ts_nd2^2 # Medida de la volatilidad. Al ser una cantidad al cuadrado, su valor ser? alto en periodos en que se experimenten grandes cambios  y comparativamente peque?o cuando sucedan cambios modestos en los precios de dichos bienes. 

plot(colones_ts_nd3)

En este gráfico de la serie ajustada por la media y elevada al cuadrado se pretende observar la volatilidad de la misma. Al ser una cantidad al cuadrado, cuando su valor ses alto en indica que se experimenten grandes cambios y comparativamente pequeño cuando sucedan cambios modestos en los precios de dichos bienes. Es así que es claro que después del 2008, alrededor del 2014 y en el 2020 es donde se presentan los mayores cambios comparativos y es efectivamente en donde se identifican crisis económicas que llevaron a la gente a buscar estas opciones de inversión.

plot(colones_ts_nd,type="l"); abline(h=0)

qqnorm(colones_ts_nd); qqline(colones_ts_nd)

acf(as.vector(colones_ts_nd))

pacf(as.vector(colones_ts_nd))

#Graficos para corroborar independencia(ruido blanco) que es diferente de correlaci?n (medida de dependencia lineal).
#Ho:residuos son independientes

acf(colones_ts_nd^2)

pacf(colones_ts_nd^2)

acf(abs(colones_ts_nd))

pacf(abs(colones_ts_nd))

En este caso algunas estacas se salen (algunas autocorrelaciones son significativas) y por tanto los rendimientos no son independientes ni identicamente distribuidos. Las autocorrelaciones significativas de los rendimientos al cuadrado o en términos absolutos reflejan la existencia de agrupamiento de volatilidad.

#McLeod.Li (Box-Ljung) test muestra una evidencia fuerte de heterocedasticidad condicional(p-value significativo). 
TSA::McLeod.Li.test(y=colones_ts_nd)

Además, a partir del test de McLeod.Li (Box-Ljung) se muestra evidencia fuerte de heterocedasticidad condicional ya que varios p-value son significativos.

3. Modelos

3.1 GARCH

garch_mod <- garchFit(~arma(0,0)+garch(1,1), data=colones_ts_nd,include.mean = FALSE)

Series Initialization:
 ARMA Model:                arma
 Formula Mean:              ~ arma(0, 0)
 GARCH Model:               garch
 Formula Variance:          ~ garch(1, 1)
 ARMA Order:                0 0
 Max ARMA Order:            0
 GARCH Order:               1 1
 Max GARCH Order:           1
 Maximum Order:             1
 Conditional Dist:          norm
 h.start:                   2
 llh.start:                 1
 Length of Series:          144
 Recursion Init:            mci
 Series Scale:              0.09442756

Parameter Initialization:
 Initial Parameters:          $params
 Limits of Transformations:   $U, $V
 Which Parameters are Fixed?  $includes
 Parameter Matrix:

 Index List of Parameters to be Optimized:
 omega alpha1  beta1 
     2      3      5 
 Persistence:                  0.9 


--- START OF TRACE ---
Selected Algorithm: nlminb 

R coded nlminb Solver: 

  0:     205.31853: 0.100000 0.100000 0.800000
  1:     205.05745: 0.113718 0.0992421 0.808603
  2:     204.71236: 0.114358 0.0841662 0.802680
  3:     204.24204: 0.137742 0.0631733 0.810655
  4:     204.18557: 0.135547 0.0478989 0.803851
  5:     203.93132: 0.151239 0.0427123 0.807213
  6:     203.89114: 0.161898 0.0318781 0.799904
  7:     203.84106: 0.173784 0.0371834 0.789182
  8:     203.10394: 0.324904 0.0990070 0.574355
  9:     202.84014: 0.469831 0.0768479 0.444926
 10:     202.73039: 0.514836 0.108307 0.371814
 11:     202.73003: 0.507096 0.119596 0.389938
 12:     202.71155: 0.498334 0.115897 0.389702
 13:     202.70982: 0.491146 0.117097 0.393875
 14:     202.70956: 0.489085 0.118707 0.394615
 15:     202.70956: 0.489077 0.118812 0.394661
 16:     202.70956: 0.489081 0.118814 0.394663

Final Estimate of the Negative LLH:
 LLH:  -137.1193    norm LLH:  -0.9522171 
      omega      alpha1       beta1 
0.004360923 0.118813794 0.394663285 

R-optimhess Difference Approximated Hessian Matrix:
             omega      alpha1       beta1
omega  -2735549.19 -15545.1878 -23464.3435
alpha1   -15545.19   -231.8720   -156.3857
beta1    -23464.34   -156.3857   -228.6424
attr(,"time")
Time difference of 0.005382061 secs

--- END OF TRACE ---


Time to Estimate Parameters:
 Time difference of 0.210856 secs

summary(garch_mod)

Title:
 GARCH Modelling 

Call:
 garchFit(formula = ~arma(0, 0) + garch(1, 1), data = colones_ts_nd, 
    include.mean = FALSE) 

Mean and Variance Equation:
 data ~ arma(0, 0) + garch(1, 1)
<environment: 0x7fed0da7baa8>
 [data = colones_ts_nd]

Conditional Distribution:
 norm 

Coefficient(s):
    omega     alpha1      beta1  
0.0043609  0.1188138  0.3946633  

Std. Errors:
 based on Hessian 

Error Analysis:
        Estimate  Std. Error  t value Pr(>|t|)  
omega   0.004361    0.001753    2.488   0.0128 *
alpha1  0.118814    0.089755    1.324   0.1856  
beta1   0.394663    0.205519    1.920   0.0548 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log Likelihood:
 137.1193    normalized:  0.9522171 

Description:
 Thu Sep 10 20:06:30 2020 by user:  


Standardised Residuals Tests:
                                Statistic p-Value     
 Jarque-Bera Test   R    Chi^2  9.301433  0.009554753 
 Shapiro-Wilk Test  R    W      0.9729107 0.00590049  
 Ljung-Box Test     R    Q(10)  33.06239  0.0002658874
 Ljung-Box Test     R    Q(15)  98.94819  2.065015e-14
 Ljung-Box Test     R    Q(20)  109.5676  2.353673e-14
 Ljung-Box Test     R^2  Q(10)  11.15519  0.3455557   
 Ljung-Box Test     R^2  Q(15)  41.28482  0.0002893284
 Ljung-Box Test     R^2  Q(20)  44.75268  0.001191913 
 LM Arch Test       R    TR^2   32.64833  0.00109813  

Information Criterion Statistics:
      AIC       BIC       SIC      HQIC 
-1.862767 -1.800896 -1.863612 -1.837627 

3.1.1 - Chequeo del modelo

acf(residuals(garch_mod)^2)

pacf(residuals(garch_mod)^2)

#1. Test de Portmanteu para residuos estandarizados al cuadrado donde la Ho es que los residuos no est?n correlacionados.

try(
gBox(garch_mod,method="absolut", plot = T)
)

Error in model$call : $ operator not defined for this S4 class

A partir de esta prueba y sumado a la inspección visual se determina que los residuos del modelo Garch estan correlacionados pues el p-value es mayor al punto de corte (0.05) y varias estacas en los gráficos son significativas.

fGarch::predict(garch_mod, n.ahead = 6,mse="uncond",plot=TRUE, crit_val=2)

NA

Por último, se observa que el ajuste del modelo en términos de predicción es muy deficiente.

3.2 Redes neuronales

nn_mod<-nnetar(colones_ts_nd)
summary(nn_mod)

          Length Class        Mode     
x         144    ts           numeric  
m           1    -none-       numeric  
p           1    -none-       numeric  
P           1    -none-       numeric  
scalex      2    -none-       list     
size        1    -none-       numeric  
subset    144    -none-       numeric  
model      20    nnetarmodels list     
nnetargs    0    -none-       list     
fitted    144    ts           numeric  
residuals 144    ts           numeric  
lags        8    -none-       numeric  
series      1    -none-       character
method      1    -none-       character
call        2    -none-       call     

acf(residuals(nn_mod)[!is.na(residuals(nn_mod))]^2)

pacf(residuals(nn_mod)[!is.na(residuals(nn_mod))]^2)

#1. Test de Portmanteu para residuos estandarizados al cuadrado donde la Ho es que los residuos no est?n correlacionados.

gBox(nn_mod,method="absolut", plot = T)

A partir de esta prueba y sumado a la inspección visual se determina que los residuos del modelo Garch estan correlacionados pues el p-value es mayor al punto de corte (0.05) y un par de estacas en los gráficos son significativas.

pred_nn_mod<-forecast::forecast(nn_mod,level = c(95), h=6, bootstrap=TRUE, npaths=10000)
pred_nn_mod

plot(pred_nn_mod)

Finalmente, se observa que aunque mejor que el modelo Garch, el modelo de redes neuronales tampoco tienen un buen ajuste en las predicciones.

3.3 Modelo autoregresivo aditivo no lineal

3.3.1 Encontrar dimensión de encrustación

dimension = estimateEmbeddingDim(colones_ts_nd, time.lag=1, max.embedding.dim=15,threshold=0.95, do.plot=TRUE)

3.3.2 Modelo

aar_mod <- aar(colones_ts_nd, m=dimension)
summary(aar_mod)

Non linear autoregressive model

AAR model

Family: gaussian 
Link function: identity 

Formula:
y ~ s(V1.0, bs = "cr") + s(V1..1, bs = "cr") + s(V1..2, bs = "cr") + 
    s(V1..3, bs = "cr") + s(V1..4, bs = "cr") + s(V1..5, bs = "cr") + 
    s(V1..6, bs = "cr") + s(V1..7, bs = "cr") + s(V1..8, bs = "cr")

Estimated degrees of freedom:
4.19 2.66 1.00 1.00 1.00 1.00 6.84 
5.27 1.00  total = 24.95 

GCV score: 0.007569904     

Residuals:
      Min        1Q    Median        3Q       Max 
-0.192080 -0.039738 -0.013724  0.038017  0.222674 

Fit:
residuals variance = 0.004716,  AIC = -607, MAPE = 1564%

Family: gaussian 
Link function: identity 

Formula:
y ~ s(V1.0, bs = "cr") + s(V1..1, bs = "cr") + s(V1..2, bs = "cr") + 
    s(V1..3, bs = "cr") + s(V1..4, bs = "cr") + s(V1..5, bs = "cr") + 
    s(V1..6, bs = "cr") + s(V1..7, bs = "cr") + s(V1..8, bs = "cr")

Parametric coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0083986  0.0067609  1.2422   0.2168

Approximate significance of smooth terms:
            edf Ref.df      F   p-value    
s(V1.0)  4.1897 4.9881 4.2915  0.001255 ** 
s(V1..1) 2.6568 3.3010 8.7961 2.236e-05 ***
s(V1..2) 1.0000 1.0000 0.8370  0.362243    
s(V1..3) 1.0000 1.0000 1.5514  0.215575    
s(V1..4) 1.0000 1.0000 0.9377  0.334995    
s(V1..5) 1.0000 1.0000 0.2848  0.594621    
s(V1..6) 6.8351 7.6297 1.0685  0.478906    
s(V1..7) 5.2686 6.0947 1.4614  0.198086    
s(V1..8) 1.0000 1.0000 2.2020  0.140692    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.301   Deviance explained = 42.6%
GCV = 0.0075699  Scale est. = 0.0061709  n = 135

plot(aar_mod)

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

NA

3.3.3 Revision del modelo

e_aar_mod <- residuals(aar_mod)
plot(e_aar_mod)

e_aar_mod <- e_aar_mod[!is.na(e_aar_mod)]
acf(e_aar_mod)

pacf(e_aar_mod)

AIC(aar_mod)

[1] -607.3787

mse(aar_mod)

[1] 0.00471599

MAPE(aar_mod)

[1] 15.64393

#fitted(aar_mod)
#coef(aar_mod)

pred_aar <- predict(aar_mod, n.ahead=6)

autoplot(ts(c(colones_ts_nd, pred_aar), start = start(colones_ts_nd), frequency = frequency(colones_ts_nd)) )

Se observa que las medidas de ajuste (MSE y MAPE) no son malas y la predicción para los 6 periodos es mejor que el modelo Garch pero hay que comparar con los demás para determinar cual obtiene mejor ajuste.

3.4 Modelo STAR (Smooth Transition AutoRegressive)

star_mod <- star(colones_ts_nd, mTh=c(0,1), control=list(maxit=10000))

Testing linearity...   p-Value =  0.1093008 
The series is linear. Using linear model instead.

summary(star_mod)

Non linear autoregressive model

AR model
Coefficients:
      const       phi.1       phi.2 
 0.01003181 -0.27505109 -0.36784039 

Residuals:
         Min           1Q       Median           3Q          Max 
-0.239698053 -0.055424235 -0.003978698  0.052357091  0.235607211 

Fit:
residuals variance = 0.007303,  AIC = -702, MAPE = 229.3%

Coefficient(s):
        Estimate  Std. Error  t value  Pr(>|t|)    
const  0.0100318   0.0072786   1.3783 0.1703347    
phi.1 -0.2750511   0.0797298  -3.4498 0.0007431 ***
phi.2 -0.3678404   0.0793699  -4.6345 8.153e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

plot(star_mod)

NA

NA

e_star_mod <- residuals(star_mod)
plot(e_star_mod)

e_star_mod <- e_star_mod[!is.na(e_star_mod)]
acf(e_star_mod)

pacf(e_star_mod)

AIC(star_mod)

[1] -702.4032

mse(star_mod)

[1] 0.007303027

MAPE(star_mod)

[1] 2.292716

pred_star <- predict(star_mod, n.ahead=6)

Los indicadores del modelo Star son mejores que los del modelo Aar, tanto el AIC, como el MSE y MAPE. Los resultados del pronóstico son similares entre los 2 modelos y es por esta razón que se procede a revisar los pronósticos de todos los modelos no lineales para determinar el mejor.

4 - Selección del mejor modelo

Una vez teniendo los modelos de suavizamiento exponencial, de regresión y Box-Jenkings (ARIMA) se procede a compararlos en términos de ajuste visual y de indicadores de ajuste para determinar cual es el mejor modelo que pronostique el número de muertes por accidentes de tránsito en Costa Rica.

# Prediccion redes neuronales
f_mod1 <-
  exp(
      log(colones_ts[145]) + cumsum(pred_nn_mod$mean)
    )

f_mod1_ts = ts(
  f_mod1,
  frequency = 12,
  start = c(2020, 2),
  end = c(2020, 7)
)

# Prediccion modelo aar
f_mod2 <-
  exp(
      log(colones_ts[145]) + cumsum(pred_aar)
    )

f_mod2_ts = ts(
  f_mod2,
  frequency = 12,
  start = c(2020, 2),
  end = c(2020, 7)
)

# Prediccion modelo Star
f_mod3 <-
  exp(
      log(colones_ts[145]) + cumsum(pred_star)
    )

f_mod3_ts = ts(
  f_mod3,
  frequency = 12,
  start = c(2020, 2),
  end = c(2020, 7)
)


todos_preds <- cbind(
  Serie = colones_ts,
  Real = colones_ts_val,
  Prediccion1 = f_mod1_ts,
  Prediccion2 = f_mod2_ts,
  Prediccion3 = f_mod3_ts
)

# Graficamos
dygraph(todos_preds, main = "Predicción todos los modelos") %>%
dySeries("Serie", label = "Cantidad") %>%
dySeries("Real", label = "Real") %>%
#dySeries("Prediccion1", label = "SES") %>%
#dySeries("Prediccion2", label = "Holt") %>%
dySeries("Prediccion1", label = "Red_neuronal") %>%
dySeries("Prediccion2", label = "AAR") %>%
dySeries("Prediccion3", label = "STAR") %>%
dyAxis("x", label = "Meses") %>% 
dyAxis("y", label = "Accidentes") %>% 
dyOptions(colors = RColorBrewer::brewer.pal(7, "Set1")) %>% 
dyRangeSelector()

acc_todos <- tibble(
  Metodo = c("Red_neuronal", "AAR", "STAR"),
  RMSE = round(
    c(
      forecast::accuracy(f_mod1_ts, colones_ts_val)[2],
      forecast::accuracy(f_mod2_ts, colones_ts_val)[2],
      forecast::accuracy(f_mod3_ts, colones_ts_val)[2]
    ),
    3
  ),
  MAE = round(
    c(
      forecast::accuracy(f_mod1_ts, colones_ts_val)[3],
      forecast::accuracy(f_mod2_ts, colones_ts_val)[3],
      forecast::accuracy(f_mod3_ts, colones_ts_val)[3]
    ),
    3
  ),
  MAPE = round(
    c(
      forecast::accuracy(f_mod1_ts, colones_ts_val)[5],
      forecast::accuracy(f_mod2_ts, colones_ts_val)[5],
      forecast::accuracy(f_mod3_ts, colones_ts_val)[5]
    ),
    3
  )
) 

acc_todos %>% 
  mutate_if(is.numeric, function(x) {
    cell_spec(x, bold = T, 
              color = spec_color(x, end = 0.9, direction = -1),
              font_size = spec_font_size(x, begin = 12,end = 14, scale_from = NA))
  }) %>%
  kable(escape = F, caption = "Medidas de ajuste (Todos los modelos)", digits = 2) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Medidas de ajuste (Todos los modelos)

Metodo	RMSE	MAE	MAPE
Red_neuronal	85345.238	66959.819	8.262
AAR	105139.802	88869.341	10.584
STAR	104337.547	88241.358	9.81

NA

Finalmente, se obtienen los valores de pronostico de 3 modelos a contrastar: Red neuronal, AAR y STAR. Se puede apreciar que el STAR fue el mas alejado de los valores reales. Entre los otros modelos se puede observar que tienen, en términos gráficos, resultados similares aunque difieren al final del periodo de prueba donde la red neuronal ajusta mejor. Moviendonos al análisis de los indicadores de ajuste se puede ver que la red neuronal fue el modelo que obtuvo los mejores indicadores y es entonces que se selecciona como el mejor modelo no lineal para predecir el número de autos importados por mes en Costa Rica.